// Copyright (C) 2014  Victor Fragoso <vfragoso@cs.ucsb.edu>
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
//     * Redistributions of source code must retain the above copyright
//       notice, this list of conditions and the following disclaimer.
//
//     * Redistributions in binary form must reproduce the above
//       copyright notice, this list of conditions and the following
//       disclaimer in the documentation and/or other materials provided
//       with the distribution.
//
//     * Neither the name of the University of California, Santa Barbara nor the
//       names of its contributors may be used to endorse or promote products
//       derived from this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL VICTOR FRAGOSO BE LIABLE FOR ANY DIRECT,
// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
// ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//

#include <vector>
#include <Eigen/Core>
#include <glog/logging.h>
#include "gtest/gtest.h"
#include "statx/distributions/evd/gpd.h"
#include "statx/distributions/evd/gpd_mle.h"

namespace libstatx {
namespace distributions {
namespace evd {

using Eigen::VectorXd;
using std::vector;

// Generated w/ Matlab: xi=1.0 sigma=1.0 mu=0.0
const vector<double> gpd_pdf {
  1, 8.264463e-01, 6.944444e-01, 5.917160e-01, 5.102041e-01, 4.444444e-01,
      3.906250e-01, 3.460208e-01, 3.086420e-01, 2.770083e-01, 2.500000e-01,
      2.267574e-01, 2.066116e-01, 1.890359e-01, 1.736111e-01, 1.600000e-01,
      1.479290e-01, 1.371742e-01, 1.275510e-01, 1.189061e-01, 1.111111e-01,
      1.040583e-01, 9.765625e-02, 9.182736e-02, 8.650519e-02, 8.163265e-02,
      7.716049e-02, 7.304602e-02, 6.925208e-02, 6.574622e-02, 6.250000e-02,
      5.948840e-02, 5.668934e-02, 5.408329e-02, 5.165289e-02, 4.938272e-02,
      4.725898e-02, 4.526935e-02, 4.340278e-02, 4.164931e-02, 4.000000e-02,
      3.844675e-02, 3.698225e-02, 3.559986e-02, 3.429355e-02, 3.305785e-02,
      3.188776e-02, 3.077870e-02, 2.972652e-02, 2.872738e-02, 2.777778e-02,
      2.687450e-02, 2.601457e-02, 2.519526e-02, 2.441406e-02, 2.366864e-02,
      2.295684e-02, 2.227668e-02, 2.162630e-02, 2.100399e-02, 2.040816e-02,
      1.983733e-02, 1.929012e-02, 1.876525e-02, 1.826150e-02, 1.777778e-02,
      1.731302e-02, 1.686625e-02, 1.643655e-02, 1.602307e-02, 1.562500e-02,
      1.524158e-02, 1.487210e-02, 1.451589e-02, 1.417234e-02, 1.384083e-02,
      1.352082e-02, 1.321178e-02, 1.291322e-02, 1.262467e-02, 1.234568e-02,
      1.207584e-02, 1.181474e-02, 1.156203e-02, 1.131734e-02, 1.108033e-02,
      1.085069e-02, 1.062812e-02, 1.041233e-02, 1.020304e-02, 1.000000e-02,
      9.802960e-03, 9.611688e-03, 9.425959e-03, 9.245562e-03, 9.070295e-03,
      8.899964e-03, 8.734387e-03, 8.573388e-03, 8.416800e-03, 8.264463e-03 };

// Generated w/ Matlab: sigma = 1.5
const vector<double> gpd_pdf_exp {
  6.666667e-01, 6.236713e-01, 5.834489e-01, 5.458205e-01, 5.106189e-01,
      4.776875e-01, 4.468800e-01, 4.180594e-01, 3.910975e-01, 3.658744e-01,
      3.422781e-01, 3.202035e-01, 2.995526e-01, 2.802336e-01, 2.621605e-01,
      2.452530e-01, 2.294359e-01, 2.146388e-01, 2.007961e-01, 1.878462e-01,
      1.757314e-01, 1.643980e-01, 1.537955e-01, 1.438767e-01, 1.345977e-01,
      1.259171e-01, 1.177963e-01, 1.101993e-01, 1.030922e-01, 9.644345e-02,
      9.022352e-02, 8.440474e-02, 7.896122e-02, 7.386877e-02, 6.910475e-02,
      6.464798e-02, 6.047864e-02, 5.657819e-02, 5.292929e-02, 4.951572e-02,
      4.632230e-02, 4.333484e-02, 4.054004e-02, 3.792549e-02, 3.547956e-02,
      3.319138e-02, 3.105077e-02, 2.904821e-02, 2.717480e-02, 2.542222e-02,
      2.378266e-02, 2.224885e-02, 2.081395e-02, 1.947160e-02, 1.821581e-02,
      1.704102e-02, 1.594200e-02, 1.491385e-02, 1.395201e-02, 1.305220e-02,
      1.221043e-02, 1.142294e-02, 1.068624e-02, 9.997051e-03, 9.352311e-03,
      8.749152e-03, 8.184893e-03, 7.657025e-03, 7.163200e-03, 6.701224e-03,
      6.269042e-03, 5.864732e-03, 5.486498e-03, 5.132657e-03, 4.801637e-03,
      4.491965e-03, 4.202264e-03, 3.931248e-03, 3.677710e-03, 3.440523e-03,
      3.218633e-03, 3.011054e-03, 2.816862e-03, 2.635194e-03, 2.465242e-03,
      2.306252e-03, 2.157514e-03, 2.018370e-03, 1.888199e-03, 1.766423e-03,
      1.652501e-03, 1.545927e-03, 1.446225e-03, 1.352954e-03, 1.265698e-03,
      1.184069e-03, 1.107705e-03, 1.036266e-03, 9.694337e-04, 9.069120e-04,
      8.484225e-04 };

// Generated w/ Matlab: xi=1.0 sigma=1.0 mu=0.0
const vector<double> gpd_cdf {
  0, 9.090909e-02, 1.666667e-01, 2.307692e-01, 2.857143e-01, 3.333333e-01,
      3.750000e-01, 4.117647e-01, 4.444444e-01, 4.736842e-01, 5.000000e-01,
      5.238095e-01, 5.454545e-01, 5.652174e-01, 5.833333e-01, 6.000000e-01,
      6.153846e-01, 6.296296e-01, 6.428571e-01, 6.551724e-01, 6.666667e-01,
      6.774194e-01, 6.875000e-01, 6.969697e-01, 7.058824e-01, 7.142857e-01,
      7.222222e-01, 7.297297e-01, 7.368421e-01, 7.435897e-01, 7.500000e-01,
      7.560976e-01, 7.619048e-01, 7.674419e-01, 7.727273e-01, 7.777778e-01,
      7.826087e-01, 7.872340e-01, 7.916667e-01, 7.959184e-01, 8.000000e-01,
      8.039216e-01, 8.076923e-01, 8.113208e-01, 8.148148e-01, 8.181818e-01,
      8.214286e-01, 8.245614e-01, 8.275862e-01, 8.305085e-01, 8.333333e-01,
      8.360656e-01, 8.387097e-01, 8.412698e-01, 8.437500e-01, 8.461538e-01,
      8.484848e-01, 8.507463e-01, 8.529412e-01, 8.550725e-01, 8.571429e-01,
      8.591549e-01, 8.611111e-01, 8.630137e-01, 8.648649e-01, 8.666667e-01,
      8.684211e-01, 8.701299e-01, 8.717949e-01, 8.734177e-01, 8.750000e-01,
      8.765432e-01, 8.780488e-01, 8.795181e-01, 8.809524e-01, 8.823529e-01,
      8.837209e-01, 8.850575e-01, 8.863636e-01, 8.876404e-01, 8.888889e-01,
      8.901099e-01, 8.913043e-01, 8.924731e-01, 8.936170e-01, 8.947368e-01,
      8.958333e-01, 8.969072e-01, 8.979592e-01, 8.989899e-01, 9.000000e-01,
      9.009901e-01, 9.019608e-01, 9.029126e-01, 9.038462e-01, 9.047619e-01,
      9.056604e-01, 9.065421e-01, 9.074074e-01, 9.082569e-01, 9.090909e-01 };

// Generated w/ Matlab: sigma = 1.5
const vector<double> gpd_cdf_exp {
  0, 6.449301e-02, 1.248267e-01, 1.812692e-01, 2.340717e-01, 2.834687e-01,
      3.296800e-01, 3.729109e-01, 4.133538e-01, 4.511884e-01, 4.865829e-01,
      5.196947e-01, 5.506710e-01, 5.796496e-01, 6.067593e-01, 6.321206e-01,
      6.558462e-01, 6.780417e-01, 6.988058e-01, 7.182307e-01, 7.364029e-01,
      7.534030e-01, 7.693068e-01, 7.841849e-01, 7.981035e-01, 8.111244e-01,
      8.233056e-01, 8.347011e-01, 8.453617e-01, 8.553348e-01, 8.646647e-01,
      8.733929e-01, 8.815582e-01, 8.891968e-01, 8.963429e-01, 9.030280e-01,
      9.092820e-01, 9.151327e-01, 9.206061e-01, 9.257264e-01, 9.305165e-01,
      9.349977e-01, 9.391899e-01, 9.431118e-01, 9.467807e-01, 9.502129e-01,
      9.534238e-01, 9.564277e-01, 9.592378e-01, 9.618667e-01, 9.643260e-01,
      9.666267e-01, 9.687791e-01, 9.707926e-01, 9.726763e-01, 9.744385e-01,
      9.760870e-01, 9.776292e-01, 9.790720e-01, 9.804217e-01, 9.816844e-01,
      9.828656e-01, 9.839706e-01, 9.850044e-01, 9.859715e-01, 9.868763e-01,
      9.877227e-01, 9.885145e-01, 9.892552e-01, 9.899482e-01, 9.905964e-01,
      9.912029e-01, 9.917703e-01, 9.923010e-01, 9.927975e-01, 9.932621e-01,
      9.936966e-01, 9.941031e-01, 9.944834e-01, 9.948392e-01, 9.951721e-01,
      9.954834e-01, 9.957747e-01, 9.960472e-01, 9.963021e-01, 9.965406e-01,
      9.967637e-01, 9.969724e-01, 9.971677e-01, 9.973504e-01, 9.975212e-01,
      9.976811e-01, 9.978307e-01, 9.979706e-01, 9.981015e-01, 9.982239e-01,
      9.983384e-01, 9.984456e-01, 9.985458e-01, 9.986396e-01, 9.987274e-01 };

// Generated w/ Matlab: xi=1.0 sigma=1.0 mu=0.0
const vector<double> gpd_quantile_data {
  0, 1.111111e-01, 2.500000e-01, 4.285714e-01, 6.666667e-01, 1,
      1.500000e+00, 2.333333e+00, 4.000000e+00, 9.000000e+00 };

// Generated w/ Matlab: xi=0.0 sigma=1.0 mu= 0.0
const vector<double> gpd_quantile_data_exp {
  0, 1.053605e-01, 2.231436e-01, 3.566749e-01, 5.108256e-01, 6.931472e-01,
      9.162907e-01, 1.203973e+00, 1.609438e+00, 2.302585e+00 };

// Rayleigh tail generated w/ Matlab: sigma=1.0 and u = 2.0
const vector<double> rayl_tail {
  1.639990e-01, 9.234254e-01, 1.209249e+00, 6.949921e-01, 3.236706e-01,
      1.250198e+00, 3.736289e-01, 7.062879e-02, 1.119897e+00, 2.007766e-01,
      4.454374e-01, 1.580042e+00, 1.788039e-01, 4.242690e-01, 3.799402e-02,
      1.455916e+00, 3.811857e-01, 1.077144e-01, 4.849757e-02, 1.524390e-01,
      1.082885e-01, 9.163793e-02, 4.149339e-01, 4.883389e-01, 4.479133e-01,
      3.473108e-02, 1.095585e-01, 1.022771e-01, 1.380078e+00, 5.192999e-01,
      1.548083e+00, 1.277406e-01, 1.074091e+00, 1.044321e-01, 6.579695e-01,
      8.099257e-02, 2.223292e-01, 6.800561e-01, 4.975733e-01, 1.224951e-01,
      1.744209e+00, 2.426436e-01, 3.884281e-03, 7.544702e-01, 4.264783e-01,
      6.507659e-01, 1.576147e-01, 6.810320e-01, 9.023873e-02, 1.105767e-02,
      3.432847e-01, 5.648951e-01, 5.772904e-01, 9.007597e-01, 2.882531e-01,
      2.003816e-01, 5.017488e-02, 5.691458e-01, 4.067266e-01, 6.390593e-01,
      8.769593e-02, 1.681791e-01, 5.371221e-01, 2.249682e-01, 2.313861e-01,
      2.217097e-01, 1.003049e+00, 4.128081e-01, 3.625543e-01, 7.606078e-01,
      1.030828e+00, 3.334478e-01, 7.718412e-01, 3.847384e-01, 2.188273e-01,
      1.057359e+00, 3.520941e-01, 2.376510e-01, 3.474343e-01, 1.230142e-01,
      2.423042e-01, 6.558300e-02, 5.360692e-02, 4.113199e-01, 2.485913e-01,
      3.343902e-01, 1.779578e+00, 1.282898e-01, 1.532888e-01, 2.781845e-01,
      9.096575e-02, 4.413923e-01, 6.116958e-01, 3.546176e-01, 5.093140e-01,
      7.313082e-03, 1.228911e-02, 6.959634e-01, 5.675502e-01, 8.988264e-02,
      3.416037e-01, 7.513501e-01, 2.294754e-01, 2.382529e-01, 9.906471e-01,
      3.688657e-01, 1.648231e-01, 8.444054e-01, 1.852471e-01, 3.297084e-01,
      1.086766e+00, 1.879264e-01, 3.818248e-01, 4.430548e-01, 4.737327e-01,
      3.275353e-01, 6.078028e-01, 1.112578e-01, 2.622672e-01, 1.085436e-01,
      7.913822e-01, 7.673774e-03, 1.029662e-01, 4.399489e-02, 6.400238e-01,
      2.080696e-01, 4.079964e-01, 1.230350e+00, 9.817702e-01, 9.034998e-01,
      1.530452e-01, 7.220132e-01, 1.321968e-01, 3.760569e-01, 4.094268e-01,
      5.586485e-01, 1.067872e-01, 2.196441e-01 };

const vector<double> rayl_tail2 {
  1.751753e-01, 1.495271e+00, 5.087389e-01, 9.188272e-01, 2.246693e-01,
      4.392266e+00, 7.225118e-01, 1.138076e+00, 4.870106e-01, 2.130642e-01,
      9.650818e-01, 2.127422e+00, 3.579166e-02, 7.035189e-01, 3.816821e-01,
      1.354670e+00, 1.655123e+00, 2.494650e+00, 1.021908e+00, 1.745268e+00,
      2.230043e+00, 2.558183e+00, 1.547089e+00, 1.296985e+00, 4.059418e-01,
      3.623963e-01, 2.465321e-01, 2.609808e-01, 5.808488e-01, 3.269070e+00,
      3.010213e+00, 1.835235e-01, 4.957769e-01, 2.137265e-02, 5.440512e-01,
      7.731908e-01, 1.106062e+00, 9.091482e-01, 6.913193e-01, 4.083195e-01,
      8.291319e-01, 3.058557e+00, 7.860067e-01, 3.141398e+00, 3.657126e+00,
      7.103165e-01, 6.374400e-02, 2.597064e-01, 6.234893e-01, 6.816143e-01,
      1.422670e+00, 3.350361e+00, 1.632757e-02, 6.095463e-01, 5.346590e-01,
      3.756810e-02, 1.161441e+00, 1.002864e+00, 9.984543e-01, 1.999213e+00,
      1.214881e+00, 5.318029e-02, 9.515359e-01, 1.602020e+00, 1.005044e+00,
      1.691282e+00, 9.099404e-01, 1.232230e+00, 2.617244e-01, 2.541753e+00,
      1.553873e+00, 1.031034e+00, 9.530544e-01, 1.913202e-01, 9.503284e-04,
      1.032037e+00, 8.373743e-01, 7.007569e-01, 1.399560e+00, 1.316709e-01,
      1.815342e+00, 3.013702e-01, 2.408118e+00, 9.190914e-01, 7.858833e-01,
      1.464256e+00, 1.175158e+00, 7.181863e-01, 1.697997e+00, 1.356154e+00,
      3.000981e+00, 1.569195e+00, 1.404160e+00, 1.424367e+00, 2.616167e+00,
      3.552943e-01, 4.317369e+00, 3.406085e+00, 4.735125e-01, 5.878241e-01,
      1.437014e-02, 4.688974e-01, 4.569186e-01, 5.075711e-01, 2.207929e+00,
      2.660913e+00, 4.432007e-01, 5.636502e-01, 2.717654e-01, 1.435360e+00,
      2.350569e+00, 2.053701e+00, 7.902972e-01, 2.364187e+00, 4.822665e-01,
      1.765282e-01, 3.446447e+00, 5.582162e-01, 6.525160e-02, 8.469489e-01,
      3.410634e-01, 4.354460e-01, 7.157257e-01, 6.967379e-01, 1.354848e+00,
      6.385138e-01, 4.116936e-01, 9.591804e-01, 1.986365e-01, 1.314003e+00,
      5.730263e-01, 6.753107e-01, 3.295620e+00, 1.431071e+00, 9.531951e-01,
      1.647478e+00, 6.541167e-01, 1.377817e+00, 2.163436e+00, 7.080685e-01,
      1.199929e-01, 1.491031e+00, 1.839319e+00, 1.060623e+00, 1.235309e-01,
      1.619875e+00, 1.732970e+00, 5.399902e-01, 1.539540e+00, 4.323152e+00,
      1.293886e-01, 5.244800e-02, 2.915442e+00, 1.855014e-01, 7.062975e-01,
      2.077772e+00, 1.308865e+00, 1.702428e+00, 6.958188e-01, 1.194509e+00,
      1.350995e+00, 2.315670e+00, 7.066359e-01, 1.639548e+00, 2.766181e+00,
      8.805990e-01, 1.745469e+00, 2.510024e+00, 2.296200e+00, 7.040282e-01,
      9.103749e-01, 1.832586e+00, 1.087482e+00, 1.143903e+00, 1.183057e+00,
      1.979991e+00, 1.865786e+00, 1.311955e+00, 2.202193e+00, 1.235208e+00,
      2.320593e+00, 2.195811e-01, 1.250008e+00, 3.021589e+00, 9.811344e-01,
      1.014605e+00, 2.713613e-01, 7.071307e-01, 1.618153e+00, 9.427606e-01,
      4.131182e-01, 9.188743e-01, 7.144681e-01, 9.280832e-01, 1.403748e+00,
      6.652099e-01, 4.499891e-01, 1.185617e+00, 1.835919e+00, 1.892322e+00,
      2.736842e+00, 1.621025e+00, 1.342167e+00, 1.171794e+00, 8.858680e-02,
      1.668492e+00, 9.145110e-01, 3.619604e+00, 7.945548e-01, 1.685852e-01,
      1.127179e-01, 1.140590e-01, 3.350129e+00, 1.115585e+00, 6.066897e+00,
      9.783046e-01, 1.646420e+00, 2.179024e-01, 1.109375e+00, 1.143569e-01,
      1.479462e+00, 5.073864e-02, 8.082524e-02, 1.298267e+00, 6.865447e-03,
      1.990075e+00, 1.669750e-01, 1.178572e+00, 3.072652e-01, 1.969641e+00,
      2.245745e+00, 1.497915e+00, 4.620586e-02, 1.411470e+00, 9.392197e-02,
      2.185391e+00, 8.889357e-01, 6.007672e-01, 2.779376e+00, 1.059026e+00,
      6.853825e-01, 3.488078e-03, 2.622136e+00, 1.682618e+00, 5.792700e-02,
      2.343236e-01, 5.663504e-01, 1.176749e+00, 1.219048e-01, 1.396612e+00,
      2.566998e-01, 2.256568e+00, 5.364467e-01, 2.446490e+00, 4.617237e-01,
      3.066177e+00, 5.883722e-01, 9.265481e-01, 1.283102e-02, 5.424989e-01,
      9.266129e-01, 6.971190e-01, 1.691642e+00, 1.184025e+00, 2.445718e+00,
      7.506538e-01, 2.755255e+00, 4.288018e+00, 2.329871e+00, 1.713017e-01 };

const vector<double> rayl_tail3 {
  2.136548e-01, 2.707218e-01, 1.119155e+00, 3.450855e+00, 1.889336e+00,
      4.897096e-01, 1.608057e-01, 9.883748e-02, 6.557958e-01, 5.929612e-01,
      1.754993e+00, 4.072330e-01, 7.401753e-01, 9.342366e-02, 3.654752e+00,
      1.647086e+00, 4.983456e-02, 1.989086e-01, 1.918696e+00, 1.698403e-01,
      3.844405e-01, 3.720623e-01, 7.439079e-01, 1.953258e+00, 5.187496e-01,
      1.682893e-02, 5.183563e-01, 1.031184e+00, 2.990411e-01, 1.704019e-01,
      1.648922e-01, 2.291928e-01, 2.535957e-01, 6.881070e-02, 2.096371e-01,
      1.638092e+00, 2.665711e-02, 1.547051e+00, 5.036355e-01, 5.702620e-01,
      2.330524e+00, 9.306992e-02, 1.733839e-01, 1.063140e+00, 1.883295e+00,
      1.073812e+00, 1.385682e+00, 2.607977e+00, 4.202758e-01, 2.820565e+00,
      1.644886e+00, 2.047088e-01, 2.427364e-01, 4.749555e-01, 1.056352e+00,
      2.668205e-01, 9.609290e-02, 9.909112e-01, 7.530018e-01, 1.055204e-01,
      2.104559e+00, 1.943403e+00, 7.329767e-01, 7.682938e-01, 1.402128e+00,
      2.009982e+00, 8.036237e-01, 4.901076e+00, 8.017449e-01, 6.604731e-01,
      1.025525e+00, 3.749865e-01, 4.311137e-01, 1.137572e+00, 6.391308e-01,
      7.198583e-01, 1.776823e+00, 1.770319e+00, 1.895467e+00, 8.461528e-02,
      1.946480e-02, 7.969864e-01, 5.526481e-01, 6.156126e-03, 1.693695e+00,
      1.872037e-01, 1.419076e+00, 1.608495e-02, 5.499211e-01, 2.178962e-01,
      2.024906e-01, 3.782849e-01, 7.317518e-03, 5.281354e-01, 1.252685e+00,
      3.987200e-01, 1.953065e-01, 6.008793e-01, 1.987154e-01, 6.323781e-01,
      5.867550e-01, 1.718075e+00, 7.685350e-01, 1.861645e+00, 1.162999e+00,
      3.375254e-01, 3.955087e+00, 7.260041e-01, 1.521279e-01, 1.523123e-01,
      3.112088e+00, 1.027333e+00, 8.931225e-01, 8.274644e-02, 1.152056e-01,
      1.175185e+00, 2.352702e-02, 1.572859e+00, 4.160540e-01, 1.953092e+00,
      1.136548e+00, 3.110009e+00, 9.787021e-01, 8.336409e-01, 3.992555e-01,
      2.275021e-01, 3.721243e-01, 4.404247e+00, 5.960064e-01, 9.605226e-01,
      1.885225e-01, 7.993633e-01, 2.007108e+00, 3.394642e+00, 4.133749e-01,
      1.118199e-02, 3.341697e-02, 1.355261e+00, 1.425273e+00, 9.031735e-01,
      6.032743e-03 };

// Generated witn xi=1.0 sigma=1.0
const vector<double> gpd_data {
  3.649519e+00, 1.312797e+00, 4.321074e-01, 1.262690e+00, 1.523320e-01,
      1.701313e+00, 3.474318e+00, 4.385702e-01, 2.443512e+00, 1.460572e-01,
      3.076213e-01, 1.088349e+00, 1.295718e+02, 3.173234e+00, 1.821959e+00,
      1.412429e-01, 5.793011e-01, 1.528987e+00, 1.750998e+00, 1.031292e+00,
      1.037962e-01, 7.938224e-02, 8.135756e-01, 5.246130e-01, 3.582283e-01,
      1.093344e+01, 6.376767e-01, 7.271343e+00, 1.076809e+00, 9.243571e-01,
      1.129557e+01, 5.840181e+00, 8.067511e-02, 7.365367e-01, 1.302307e+00,
      1.315213e+00, 3.140494e+00, 3.170728e+00, 1.376750e+00, 1.450502e+00,
      2.587296e-01, 8.666972e-01, 9.335795e-01, 3.520471e+00, 1.243067e-01,
      4.351209e-01, 1.884934e-01, 4.684062e-01, 1.627966e-02, 2.052172e-02,
      1.998776e+00, 5.306914e-01, 3.436460e-01, 8.336165e-01, 1.695237e-01,
      2.577268e+00, 3.676328e+00, 3.586126e+01, 3.113889e+00, 7.364907e-01,
      7.988221e-01, 1.222389e+00, 1.588655e+00, 3.154153e+01, 2.150657e-01,
      5.030947e-01, 9.869404e-01, 2.048170e+00, 2.313521e-01, 2.033022e-01,
      1.665146e-01, 9.258666e-01, 2.996108e-02, 1.805648e-01, 5.136115e+00,
      2.103059e+00, 2.119092e-01, 1.139486e+00, 9.288127e-01, 5.029963e-02,
      3.594080e-02, 2.862222e-01, 6.385923e-01, 8.977742e-02, 8.414354e-01,
      2.908624e+00, 1.502860e+00, 1.673368e+00, 6.766819e-01, 5.699327e+00,
      5.615092e+00, 2.897247e+00, 5.472893e-02, 2.207307e-01, 3.728020e-01,
      3.043534e-02, 1.575069e-03, 1.266448e-02, 3.162488e+00, 4.286063e-02 };

TEST(GPD, PDF_Exponential) {
  // Exponential case
  const double xi = 0.0;
  const double sigma = 1.5;
  const double dx = 0.1;
  double x = 0.0;
  for (int i = 0; i < gpd_pdf_exp.size(); i++) {
    double y_gt = gpd_pdf_exp[i];
    double y = gpdpdf(x, xi, sigma);
    EXPECT_NEAR(y_gt, y, 1e-3);
    x += dx;
  }
}

TEST(GPD, PDF_NonExponential) {
  // NonExponential case
  const double xi = 1.0;
  const double sigma = 1.0;
  const double dx = 0.1;
  double x = 0.0;
  for (int i = 0; i < gpd_pdf.size(); i++) {
    double y_gt = gpd_pdf[i];
    double y = gpdpdf(x, xi, sigma);
    EXPECT_NEAR(y_gt, y, 1e-3);
    x += dx;
  }
}

TEST(GPD, CDF_Exponential) {
  // Exponential case
  const double xi = 0.0;
  const double sigma = 1.5;
  const double dx = 0.1;
  double x = 0.0;
  for (int i = 0; i < gpd_cdf_exp.size(); i++) {
    double y_gt = gpd_cdf_exp[i];
    double y = gpdcdf(x, xi, sigma);
    EXPECT_NEAR(y_gt, y, 1e-3);
    x += dx;
  }
}

TEST(GPD, CDF_NonExponential) {
  // Non-Exponential case
  const double xi = 1.0;
  const double sigma = 1.0;
  const double dx = 0.1;
  double x = 0.0;
  for (int i = 0; i < gpd_cdf.size(); i++) {
    double y_gt = gpd_cdf[i];
    double y = gpdcdf(x, xi, sigma);
    EXPECT_NEAR(y_gt, y, 1e-3);
    x += dx;
  }
}

TEST(GPD, InvalidParams) {
  // Non exponential -- out of domain
  EXPECT_EQ(std::numeric_limits<double>::infinity(), gpdpdf(50.0, 1.0, -1.0));
  // Invalid sigma
  EXPECT_EQ(std::numeric_limits<double>::infinity(), gpdpdf(1.0, 1.0, 0.0));
}

TEST(GPD, FitNoData) {
  vector<double> data;
  double xi, sigma;
  EXPECT_FALSE(gpdfit(data, &xi, &sigma));
}

TEST(GPD, Quantile_NonExponential) {
  const double sigma = 1.0;
  const double xi = 1.0;
  double p = 0.0;
  const double dp = 0.1;
  for (int i = 0; i < gpd_quantile_data.size(); i++) {
    double q_gt = gpd_quantile_data[i];
    double q = gpd_quantile(p, xi, sigma);
    EXPECT_NEAR(q_gt, q, 1e-3);
    p += dp;
  }
}

TEST(GPD, Quantile_Exponential) {
  const double sigma = 1.0;
  const double xi = 0.0;
  double p = 0.0;
  const double dp = 0.1;
  for (int i = 0; i < gpd_quantile_data_exp.size(); i++) {
    double q_gt = gpd_quantile_data_exp[i];
    double q = gpd_quantile(p, xi, sigma);
    EXPECT_NEAR(q_gt, q, 1e-3);
    p += dp;
  }
}

TEST(GPD, MLE_Objective) {
  const double mle_gt = 25.0869;
  GPDMLEObjective mle(rayl_tail);
  // GPD params: xi=-0.1730 sigma=0.5246
  VectorXd x(2);
  x(0) = -0.1730;  // xi
  x(1) = 0.5246;  // sigma
  double mle_val = mle(x);
  EXPECT_NEAR(mle_gt, mle_val, 1e-3);
}


TEST(GPD, MLE_Objective2) {
  // GPD params: xi=-0.1870  sigma=1.4402
  const double mle_gt = 318.0067;
  GPDMLEObjective mle(rayl_tail2);
  VectorXd x(2);
  x(0) = -0.1870;  // xi
  x(1) = 1.4402;  // sigma
  double mle_val = mle(x);
  EXPECT_NEAR(mle_gt, mle_val, 1e-3);
}

TEST(GPD, MLE_Objective3) {
  // GPD params: xi=0.0247 sigma=0.9108
  const double mle_gt = 131.3145;
  GPDMLEObjective mle(rayl_tail3);
  VectorXd x(2);
  x(0) = 0.0247;  // xi
  x(1) = 0.9108;  // sigma
  double mle_val = mle(x);
  EXPECT_NEAR(mle_gt, mle_val, 1e-3);
}

TEST(GPD, MLE_Gradient1) {
  const int n = 2;
  // Testcase 1
  GPDMLEGradientFunctor gradient(rayl_tail);
  VectorXd x(n), g(n);
  x(0) = -0.1730;  // xi
  x(1) = 0.5246;  // sigma
  gradient(x, &g);
  VLOG(1) << "Gradient: " << g.transpose();
  EXPECT_LT(g.norm(), sqrt(n));
}

TEST(GPD, MLE_Gradient2) {
  const int n = 2;
  VectorXd x(n), g(n);
  // Testcase 2
  GPDMLEGradientFunctor gradient(rayl_tail2);
  x(0) = -0.1870;  // xi
  x(1) = 1.4402;  // sigma
  gradient(x, &g);
  VLOG(1) << "Gradient: " << g.transpose();
  EXPECT_LT(g.norm(), sqrt(n));
}

TEST(GPD, MLE_Gradient3) {
  const int n = 2;
  VectorXd x(n), g(n);
  // Testcase 3
  GPDMLEGradientFunctor gradient(rayl_tail3);
  x(0) = 0.0247;  // xi
  x(1) = 0.9108;  // sigma
  gradient(x, &g);
  VLOG(1) << "Gradient: " << g.transpose();
  EXPECT_LT(g.norm(), sqrt(n));
}

TEST(GPD, FitMLE1) {
  double xi, sigma;  // Unknowns
  const int n = 2;
  // Case 1
  VectorXd x(n);
  x(0) = -0.1730;  // xi
  x(1) = 0.5246;  // sigma
  VLOG(1) << "Ground Truth: " << x.transpose();
  EXPECT_TRUE(gpdfit(rayl_tail, &xi, &sigma));
  VLOG(1) << "sigma: " << sigma << " xi=" << xi;
  EXPECT_NEAR(x(0), xi, 1.0);
  EXPECT_NEAR(x(1), sigma, 1.0);
}

TEST(GPD, FitMLE2_BadInitialPoint) {
  double xi, sigma;  // Unknowns
  const int n = 2;
  VectorXd x(n);
  x(0) = -0.1870;  // xi
  x(1) = 1.4402;  // sigma
  VLOG(1) << "Ground Truth: " << x.transpose();
  EXPECT_FALSE(gpdfit(rayl_tail2, &xi, &sigma));
  VLOG(1) << "sigma: " << sigma << " xi=" << xi;
}

TEST(GPD, FitMLE3) {
  double xi, sigma;  // Unknowns
  const int n = 2;
  VectorXd x(n);
  x(0) = 0.0247;  // xi
  x(1) = 0.9108;  // sigma
  VLOG(1) << "Ground Truth: " << x.transpose();
  EXPECT_TRUE(gpdfit(rayl_tail3, &xi, &sigma));
  VLOG(1) << "sigma: " << sigma << " xi=" << xi;
  EXPECT_NEAR(x(0), xi, 1.0);
  EXPECT_NEAR(x(1), sigma, 1.0);
}

TEST(GPD, FitMLE4) {
  // TEST with R!
  double xi, sigma;  // Unknowns
  const int n = 2;
  VectorXd x(n);
  x(0) = 0.5;  // xi
  x(1) = 1.0;  // sigma
  VLOG(1) << "Ground Truth: " << x.transpose();
  EXPECT_TRUE(gpdfit(gpd_data, &xi, &sigma));
  VLOG(1) << "Computed: " << xi << " " << sigma;
  EXPECT_NEAR(x(0), xi, 1.0);
  EXPECT_NEAR(x(1), sigma, 1.0);
}
#ifdef STATX_WITH_CERES
//////////////////////////////////////////
// Testing GEV fit w/ CERES
TEST(GPD, FitQuantileLeastSquares_Case1) {
  double xi, sigma;  // Unknowns
  const int n = 2;
  // Case 1
  VectorXd x(n);
  x(0) = -0.1730;  // xi
  x(1) = 0.5246;  // sigma
  VLOG(1) << "Ground Truth: " << x.transpose();
  EXPECT_TRUE(gpdfit(rayl_tail, &xi, &sigma, QUANTILE_NLS));
  VLOG(1) << "sigma: " << sigma << " xi=" << xi;
  EXPECT_NEAR(x(0), xi, 1.0);
  EXPECT_NEAR(x(1), sigma, 1.0);
}

TEST(GPD, FitQuantileLeastSquares_Case2) {
  // TEST with R!
  double xi, sigma;  // Unknowns
  const int n = 2;
  VectorXd x(n);
  x(0) = 0.5;  // xi
  x(1) = 1.0;  // sigma
  VLOG(1) << "Ground Truth: " << x.transpose();
  EXPECT_TRUE(gpdfit(gpd_data, &xi, &sigma, QUANTILE_NLS));
  VLOG(1) << "Computed: " << xi << " " << sigma;
  EXPECT_NEAR(x(0), xi, 1.0);
  EXPECT_NEAR(x(1), sigma, 1.0);
}

TEST(GPD, FitQuantileLeastSquares_Case3) {
  double xi, sigma;  // Unknowns
  const int n = 2;
  VectorXd x(n);
  x(0) = 0.0247;  // xi
  x(1) = 0.9108;  // sigma
  VLOG(1) << "Ground Truth: " << x.transpose();
  EXPECT_TRUE(gpdfit(rayl_tail3, &xi, &sigma, QUANTILE_NLS));
  VLOG(1) << "sigma: " << sigma << " xi=" << xi;
  EXPECT_NEAR(x(0), xi, 1.0);
  EXPECT_NEAR(x(1), sigma, 1.0);
}
#endif
}  // evd
}  // distributions
}  // statx
